5: Compute Standardize Evotranspiration Index based on spatial point


1 Introduction

In this tutorial, we extract the CHIRPS precipitation observations in Suriname, based on the centroids of the Principal Sampling Units (PSU) of the Suriname Survey of Living Conditions for the years 2016/17 and 2022. Based on the extracted precipitation, the tutorial shows how to compute the Standardized Precipitation Index (SPI) and merge it with the survey based on the date of interviews.

2 Code

2.1 Set up

We start by setting up the stage for our analysis.

First, we load the necessary packages. We load only climatic4economist package that contains several functions meant to extract and merge spatial variables with surveys. During the tutorial we will use other packages but instead of loading all the package at the begging we will call specific function each time.

library(climatic4economist)

In the setup, we also want to create the paths to the various data sources and load the necessary functions for extraction. Note .. means one step back to the folder directory, i.e. one folder back.

Note that how to set up the paths depends on your folder organization but there are overall two approaches:

  1. you can use the R project, by opening the project directly you don’t need to set up the path to the project. Automatically the project figures out on its own where it is located in the computer and set that path as working folder.
  2. you can manually set the working folder with the function setwd().
# path to data folder
path_to_data <- file.path("..",
                          "..", "data")

# survey 
path_to_survey <- file.path(path_to_data, "survey", "AFR")

# administrative division
path_to_adm_div <- file.path(path_to_data, "adm_div", "geoBoundaries")

# weather variables
path_to_pre_monthly <- file.path(path_to_data, "weather", "ERA5_Land", "monthly_pre",
                                 "era5_monthly_avg_total_pre_50_25.nc")
path_to_evo_monthly <- file.path(path_to_data, "weather", "ERA5_Land", "monthly_evo",
                                 "era5_monthly_total_evo_50_25.nc")
# to result folder
path_to_result <- file.path(path_to_data, "result")
1
concatenate the string to make a path
2
.. means one folder back


2.2 Read the data

2.2.1 Survey

We start by reading the surveys data. The surveys must have an id that uniquely identifies the household and the coordinates of their interviews.

The next steps are a bit convoluted. Lets split one by one.

  1. list.files(path_to_survey, full.names = TRUE). The surveys are stored in two files. We use list.files() to list the files.
  2. lapply(haven::read_dta). We use the lapply() to apply to each of them the function haven::read_dta(). This last function actually read dta files into R. The result are two separate block of data, each for each file. They are two separated item within a list.
  3. At this point we can bind the two separate data into a single one with the function dplyr::bind_rows(). We could have done it before but like this we ensured the two data sets have the same columns.
srvy <- list.files(path_to_survey, pattern = "dta$", full.names = TRUE) |>
    lapply(haven::read_dta) |>
    dplyr::bind_rows() |> 
    dplyr::filter(country != "Uganda") |>
    # dplyr::select(dplyr::any_of(c("country",  "year", "hhid", "lat", "lon", "interview_date"))) |>
    dplyr::group_by(country) |>
    dplyr::mutate(interview_date = clock::date_parse(interview_date, format = "%d/%m/%Y"),
                  survey_year = clock::get_year(median(interview_date, na.rm = TRUE)),
                  survey = paste0(country, substr(survey_year, 3, 4)),
                  survey = gsub(" ", "", survey),
                  .before = hhid)

surveys <- srvy$survey |> unique() |> sort()
surveys
[1] "BurkinaFaso22" "Ethiopia19"    "Malawi19"      "Mali22"       
[5] "Nigeria18"     "Tanzania21"    "Togo22"       

We can keep all the surveys together but the size can be a challenge. Therefore, we split them with the function dplyr::group_split(). Now, each survey is a separate block, but they are all stored in the same list. We also distinguish the surveys that have spatial coordinates from the others.

no_coord_srvy <- c("Nigeria18", "Tanzania21")
srvy_adm_cntry <- srvy |>
    dplyr::filter(survey %in% no_coord_srvy) |> 
    dplyr::select(-c(ea_id, TA_code, ea_code, 
                     zone_code, city_code, subcity_code)) |>
    dplyr::mutate(dplyr::across(.cols = dplyr::where(is.character),
                                .fns = ~ifelse(.x == "", NA_character_, .x)),
                  dplyr::across(.cols = dplyr::where(labelled::is.labelled),
                                .fns = labelled::to_character)) |>
    dplyr::group_by(survey) |>
    dplyr::group_split() |>
    purrr::map(janitor::remove_empty,
               which = "cols") |>
    setNames(no_coord_srvy) |> 
  purrr::map(\(x)  x |> 
               dplyr::mutate(dplyr::across(.cols = dplyr::where(is.character),
                                           .fns = ~ ifelse(.x == "" | grepl("CONFI", .x),
                                                           NA_character_, .x)) )) |> 
  purrr::map(janitor::remove_empty, which = "cols")

srvy_adm_cntry
$Nigeria18
# A tibble: 40,183 × 17
   survey_year hhid  country admin1     admin2 admin3   lat   lon interview_date
         <int> <chr> <chr>   <chr>      <chr>  <chr>  <dbl> <dbl> <date>        
 1        2018 10001 Nigeria South East Abia   Abia      NA    NA 2018-08-20    
 2        2018 10002 Nigeria South East Abia   Abia      NA    NA 2018-08-21    
 3        2018 10003 Nigeria South East Abia   Abia      NA    NA 2018-08-20    
 4        2018 10004 Nigeria South East Abia   Abia      NA    NA 2018-08-21    
 5        2018 10005 Nigeria South East Abia   Abia      NA    NA 2018-08-21    
 6        2018 10008 Nigeria South East Abia   Abia      NA    NA 2018-08-21    
 7        2018 10009 Nigeria South East Abia   Abia      NA    NA 2018-08-22    
 8        2018 10010 Nigeria South East Abia   Abia      NA    NA 2018-08-20    
 9        2018 10031 Nigeria South East Abia   Abia      NA    NA 2018-08-07    
10        2018 10032 Nigeria South East Abia   Abia      NA    NA 2018-08-08    
# ℹ 40,173 more rows
# ℹ 8 more variables: survey <chr>, admin1_code <chr>, admin2_code <chr>,
#   admin3_code <chr>, admin1_desc <chr>, admin2_desc <chr>, admin3_desc <chr>,
#   year <chr>

$Tanzania21
# A tibble: 4,709 × 11
   survey_year hhid        country interview_date survey region_name region_code
         <int> <chr>       <chr>   <date>         <chr>  <chr>             <dbl>
 1        2021 1000-001-01 Tanzan… 2021-09-15     Tanza… arusha                2
 2        2021 1000-001-02 Tanzan… 2021-09-15     Tanza… arusha                2
 3        2021 1000-001-03 Tanzan… 2021-09-15     Tanza… arusha                2
 4        2021 1000-001-06 Tanzan… 2021-09-15     Tanza… arusha                2
 5        2021 1001-001-01 Tanzan… 2021-09-15     Tanza… arusha                2
 6        2021 1002-001-01 Tanzan… 2021-09-15     Tanza… arusha                2
 7        2021 1003-001-01 Tanzan… 2021-09-15     Tanza… arusha                2
 8        2021 1005-001-01 Tanzan… 2021-09-15     Tanza… arusha                2
 9        2021 1006-001-01 Tanzan… 2021-09-15     Tanza… arusha                2
10        2021 1006-001-03 Tanzan… 2021-09-15     Tanza… arusha                2
# ℹ 4,699 more rows
# ℹ 4 more variables: district_name <chr>, district_code <dbl>,
#   ward_code <dbl>, village_code <dbl>
ys_coord_srvy <- setdiff(surveys, no_coord_srvy)

srvy_coord_cntry <- srvy |>
    dplyr::filter(survey != "Nigeria18") |>
    dplyr::filter(!is.na(lat) & lat > -999999999) |>
    dplyr::filter(!(lat == 0.00000 & lon == 0.00000)) |>
    dplyr::select(c(survey, country, hhid, lat, lon, interview_date)) |>
    dplyr::group_by(survey) |>
    dplyr::group_split() |>
    setNames(ys_coord_srvy)

srvy_coord_cntry
<list_of<
  tbl_df<
    survey        : character
    country       : character
    hhid          : character
    lat           : double
    lon           : double
    interview_date: date
  >
>[5]>
$BurkinaFaso22
# A tibble: 3,227 × 6
   survey        country      hhid     lat     lon interview_date
   <chr>         <chr>        <chr>  <dbl>   <dbl> <date>        
 1 BurkinaFaso22 Burkina Faso 002001  14.4 -0.231  2021-12-21    
 2 BurkinaFaso22 Burkina Faso 002002  14.4 -0.231  2021-12-21    
 3 BurkinaFaso22 Burkina Faso 002004  14.4 -0.231  2021-12-21    
 4 BurkinaFaso22 Burkina Faso 002005  14.4 -0.231  2021-12-21    
 5 BurkinaFaso22 Burkina Faso 002007  14.4 -0.231  2021-12-21    
 6 BurkinaFaso22 Burkina Faso 002040  14.5 -0.234  2021-12-22    
 7 BurkinaFaso22 Burkina Faso 003008  14.0 -0.0260 2022-06-29    
 8 BurkinaFaso22 Burkina Faso 003019  12.4 -1.54   2022-07-25    
 9 BurkinaFaso22 Burkina Faso 003021  14.0 -0.0320 2022-06-26    
10 BurkinaFaso22 Burkina Faso 003027  14.0 -0.0260 2022-06-22    
# ℹ 3,217 more rows

$Ethiopia19
# A tibble: 17,830 × 6
   survey     country  hhid                 lat   lon interview_date
   <chr>      <chr>    <chr>              <dbl> <dbl> <date>        
 1 Ethiopia19 Ethiopia 010101088800910007  14.3  37.8 2019-06-19    
 2 Ethiopia19 Ethiopia 010101088800910017  14.3  37.8 2019-06-18    
 3 Ethiopia19 Ethiopia 010101088800910026  14.3  37.8 2019-06-18    
 4 Ethiopia19 Ethiopia 010101088800910029  14.3  37.8 2019-06-17    
 5 Ethiopia19 Ethiopia 010101088800910038  14.3  37.8 2019-06-07    
 6 Ethiopia19 Ethiopia 010101088800910046  14.3  37.8 2019-06-04    
 7 Ethiopia19 Ethiopia 010101088800910054  14.3  37.8 2019-06-05    
 8 Ethiopia19 Ethiopia 010101088800910062  14.3  37.8 2019-06-08    
 9 Ethiopia19 Ethiopia 010101088800910070  14.3  37.8 2019-06-17    
10 Ethiopia19 Ethiopia 010101088800910082  14.3  37.8 2019-06-06    
# ℹ 17,820 more rows

$Malawi19
# A tibble: 16,955 × 6
   survey   country hhid           lat   lon interview_date
   <chr>    <chr>   <chr>        <dbl> <dbl> <date>        
 1 Malawi19 Malawi  101011000014 -9.70  33.2 2019-08-29    
 2 Malawi19 Malawi  101011000023 -9.70  33.2 2019-08-29    
 3 Malawi19 Malawi  101011000040 -9.70  33.2 2019-08-28    
 4 Malawi19 Malawi  101011000071 -9.70  33.2 2019-08-29    
 5 Malawi19 Malawi  101011000095 -9.70  33.2 2019-08-28    
 6 Malawi19 Malawi  101011000115 -9.70  33.2 2019-08-29    
 7 Malawi19 Malawi  101011000126 -9.70  33.2 2019-08-28    
 8 Malawi19 Malawi  101011000135 -9.70  33.2 2019-08-29    
 9 Malawi19 Malawi  101011000183 -9.70  33.2 2019-08-28    
10 Malawi19 Malawi  101011000190 -9.70  33.2 2019-08-28    
# ℹ 16,945 more rows

$Mali22
# A tibble: 6,119 × 6
   survey country hhid    lat   lon interview_date
   <chr>  <chr>   <chr> <dbl> <dbl> <date>        
 1 Mali22 Mali    00101  14.7 -12.2 2021-11-10    
 2 Mali22 Mali    00102  14.7 -12.2 2021-11-10    
 3 Mali22 Mali    00103  14.7 -12.2 2021-11-11    
 4 Mali22 Mali    00104  14.7 -12.2 2021-11-11    
 5 Mali22 Mali    00105  14.7 -12.2 2021-11-10    
 6 Mali22 Mali    00106  14.7 -12.2 2021-11-10    
 7 Mali22 Mali    00107  14.7 -12.2 2021-11-10    
 8 Mali22 Mali    00108  14.7 -12.2 2021-11-10    
 9 Mali22 Mali    00109  14.7 -12.2 2021-11-10    
10 Mali22 Mali    00110  14.7 -12.2 2021-11-12    
# ℹ 6,109 more rows

$Togo22
# A tibble: 6,462 × 6
   survey country hhid     lat   lon interview_date
   <chr>  <chr>   <chr>  <dbl> <dbl> <date>        
 1 Togo22 Togo    001002  6.16  1.23 2021-11-13    
 2 Togo22 Togo    001004  6.16  1.23 2021-11-14    
 3 Togo22 Togo    001005  6.16  1.23 2021-11-15    
 4 Togo22 Togo    001006  6.16  1.23 2021-11-15    
 5 Togo22 Togo    001007  6.16  1.23 2021-11-14    
 6 Togo22 Togo    001008  6.16  1.23 2021-11-13    
 7 Togo22 Togo    001009  6.16  1.23 2021-11-14    
 8 Togo22 Togo    001010  6.16  1.23 2021-11-15    
 9 Togo22 Togo    001011  6.16  1.23 2021-11-15    
10 Togo22 Togo    001013  6.16  1.23 2021-11-15    
# ℹ 6,452 more rows

2.2.2 Weather data

To weather data is stored as nc files. We can read them with the function terra::rast().

Note how all the data sets have the same coordinate reference system (CRS), i.e. EPSG:4326. This is important because in this way all the data can “spatially” talk to each other.

tpre <- terra::rast(path_to_pre_monthly)
names(tpre) <- terra::names(tpre) |> second_to_date()

tevo <- terra::rast(path_to_evo_monthly)
names(tevo) <- terra::names(tevo) |> second_to_date()

2.2.3 Administrative Boundaries

We now move to read the administrative divisions. We use the function read_adm_div() to do so. This function looks for spatial polygons for the iso and lvl provided provided.

Even if we have the coordinates from the survey, we will extract some spatial variables at the administrative division.

adm_div <- read_geoBoundaries(path_to_adm_div,
                              iso = c("BFA", "ETH", "MLI", "MWI", "TZA", "TGO"),
                              lvl = 2) |>
  setNames(c("BurkinaFaso22", "Ethiopia19", "Mali22", "Malawi19",
             "Togo22", "Tanzania21"))

Nigeria <- read_geoBoundaries(path_to_adm_div,
                              iso = c("NGA"),
                              lvl = 1)
adm_div$Nigeria18 <- Nigeria

adm_div <- adm_div[sort(names(adm_div))]

adm_div$BurkinaFaso22
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 45, 4  (geometries, attributes)
 extent      : -5.51892, 2.4054, 9.40111, 15.08259  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       : ID_adm_div   iso         adm_div_1 adm_div_2
 type        :      <chr> <chr>             <chr>     <chr>
 values      :          1   BFA Boucle du Mouhoun      Bale
                        2   BFA Boucle du Mouhoun   Mouhoun
                        3   BFA Boucle du Mouhoun    Nayala

2.3 Georeference the Surveys

srvy_coord_geo <- srvy_coord_cntry |> 
  list(purrr::keep_at(adm_div, ys_coord_srvy), names(srvy_coord_cntry)) |>
  purrr::pmap(get_poly_attr_for_point)
srvy_adm_geo <- srvy_adm_cntry
srvy_adm_geo$Nigeria18 <- srvy_adm_cntry$Nigeria18 |> 
  dplyr::mutate(adm_div_1 = stringr::str_to_title(admin1),
                adm_div_1 = dplyr::case_when(
                  adm_div_1 == "Fct" ~ "Abuja Federal Capital Territory",
                  .default = adm_div_1))

srvy_adm_geo$Tanzania21 <- srvy_adm_cntry$Tanzania21 |> 
  dplyr::mutate(adm_div_1 = stringr::str_to_title(region_name),
                adm_div_2 = stringr::str_to_title(district_name),
                
                adm_div_1 = dplyr::case_when(
                  adm_div_1 == "Dar Es Salaam" ~ "Dar es Salaam",
                  adm_div_1 == "Kaskazini Pemba" ~ "North Pemba",
                  adm_div_1 == "Kaskazini Unguja" ~ "Zanzibar North",
                  adm_div_1 == "Kusini Pemba" ~ "South Pemba",
                  adm_div_1 == "Kusini Unguja" ~ "Zanzibar South & Central",
                  adm_div_1 == "Mjini Magharibi Unguja" ~ "Zanzibar Urban/West",
                  # Songwe was created the 2016 from the western half of Mbeya Region
                  adm_div_1 == "Songwe" ~ "Mbeya", 
                  
                  # wrong adm_div based on district
                  adm_div_2 == "Masasi Rural" ~ "Mtwara", 
                  adm_div_2 == "Kigamboni" ~ "Dar es Salaam", 
                  grepl("Kahama", adm_div_2) & adm_div_1 == "Katavi" ~ "Shinyanga",  # !!!
                  adm_div_2 == "Chakechake" ~ "South Pemba", 
                  adm_div_2 == "Kinondoni" ~ "Dar es Salaam", 
                  adm_div_2 == "Songea" ~ "Ruvuma",
                  grepl("Songea", adm_div_2) ~ "Ruvuma",
              
                  adm_div_2 == "Ngorongoro" ~ "Arusha", 
                  adm_div_2 == "Wete" ~ "North Pemba", 
                  adm_div_2 == "Tandahimba" ~ "Mtwara", 
                  adm_div_2 == "Babati" ~ "Manyara",
                  grepl("Baba", adm_div_2) ~ "Manyara", # !!!!
                  grepl("Nzega", adm_div_2) ~ "Tabora",
                  adm_div_2 == "Makete" ~ "Njombe",
                  adm_div_2 == "Mbeya Urban" ~ "Mbeya",
                  .default = adm_div_1),
                
                adm_div_2 = dplyr::case_when(
                  # adm_div_1  Arusha
                  adm_div_2 == "Arusha Rural" ~ "Arusha",
                  # adm_div_1  Shinyanga
                  adm_div_2 == "Kahama Rural" ~ "Kahama",
                  adm_div_2 == "Kahama Town" ~ "Kahama Township Authority",
                  adm_div_2 == "Shinyanga Rural" ~ "Shinyanga",
                  adm_div_2 == "Nzega Town" ~ "Nzega",
                  # in 2012 by splitting the Kahama District into Msalala and Ushetu
                  adm_div_2 == "Msalala" ~ "Kahama",
                  adm_div_2 == "Ushetu" ~ "Kahama",
                  # adm_div_1  Katavi
                  adm_div_2 == "Mpanda Rural" ~ "Mpanda",
                  # established in 2012 from Mlele
                  adm_div_2 == "Mpimbwe" ~ "Mlele",
                  adm_div_2 == "Nsimbo" ~ "Mlele",
                  # adm_div_1  Singida
                  adm_div_2 == "Singida Rural" ~ "Singida",
                  # established in 2015 from Manyoni
                  adm_div_2 == "Itigi" ~ "Manyoni",
                  # adm_div_1 Dar es Salaam
                  # In 2015 Temeke was divided into Temeke and Kigamboni
                  adm_div_2 == "Kigamboni" ~ "Temeke",
                  # Kinondoni should be onlyu the eastern part
                  adm_div_2 == "Ubungo" ~ "Kinondoni",
                  # adm_div_1 Kagera
                  adm_div_2 == "Bukoba Rural" ~ "Bukoba",
                  # adm_div_1 Dodoma
                  adm_div_2 == "Kondoa Urban" ~ "Kondoa",
                  # adm_div_1 Mbeya
                  adm_div_2 == "Mbeya Rural" ~ "Mbeya",
                  adm_div_2 == "Mbalali" ~ "Mbarali",
                  # created in 2013, not perfect geo match
                  adm_div_2 == "Busokelo" ~ "Makete",
                  # adm_div_1 Morogoro
                  adm_div_2 == "Ifakara Urban" ~ "Kilombero",
                  adm_div_2 == "Morogoro Rural" ~ "Morogoro",
                  # it should be the western part of Ulanga
                  adm_div_2 == "Malinyi" ~ "Ulanga",
                  # adm_div_1 Kigoma
                  adm_div_2 == "Kasulu Rural" ~ "Kasulu",
                  adm_div_2 == "Kasulu Town" ~ "Kasulu Township Authority",
                  adm_div_2 == "Kigoma Rural" ~ "Kigoma",
                  adm_div_2 == "Kigoma Ujiji Urban" ~ "Kigoma  Urban",
                  # adm_div_1 Mtwara
                  adm_div_2 == "Masasi Rural" ~ "Masasi",
                  adm_div_2 == "Masasi Urban" ~ "Masasi  Township Authority",
                  adm_div_2 == "Mtwara Mikindani" ~ "Mtwara Urban",
                  adm_div_2 == "Mtwara Rural" ~ "Mtwara",
                  # adm_div_1 Geita
                  adm_div_2 == "Geita Town" ~ "Geita",
                  # adm_div_1 Mwanza
                  adm_div_2 == "Mwanza Urban" ~ "Nyamagana",
                  # created in 2015, from the eastern part of Sengerema
                  adm_div_2 == "Buchosa" ~ "Sengerema",
                  # adm_div_1 Iringa
                  adm_div_2 == "Iringa Rural" ~ "Iringa",
                  adm_div_2 == "Mafinga Town" ~ "Mafinga Township Authority",
                  # adm_div_1 Njombe
                  adm_div_2 == "Makambako Town" ~ "Makambako Township Authority",
                  adm_div_2 == "Mbeya Urban" ~ "Mbeya",
                  adm_div_2 == "Njombe Rural" ~ "Njombe",
                  adm_div_2 == "Njombe Town" ~ "Njombe Urban",
                  # adm_div_1 South Pemba
                  adm_div_2 == "Chakechake" ~ "Chake Chake",
                  # adm_div_1 Lindi
                  adm_div_2 == "Lindi Rural" ~ "Lindi",
                  # adm_div_1 Manyara
                  adm_div_2 == "Babati Rural" ~ "Babati",
                  adm_div_2 == "Babati Town" ~ "Babati UrbanBabati Urban",
                  adm_div_2 == "Mbulu Town" ~ "Mbulu",
                  # adm_div_1 Mara
                  adm_div_2 == "Butiama" ~ "Babati",
                  adm_div_2 == "Musoma Rural" ~ "Musoma",
                  # adm_div_1 Ruvuma
                  adm_div_2 == "Songea Rural" ~ "Songea",
                  # not clear
                  adm_div_2 == "Madaba" ~ "Songea",
                  # adm_div_1 Simiyu
                  adm_div_2 == "Bariadi Town" ~ "Bariadi",
                  # adm_div_1 Rukwa
                  adm_div_2 == "Sumbawanga Rural" ~ "Sumbawanga",
                  # adm_div_1 Pwani
                  adm_div_2 == "Kibaha Rural" ~ "Kibaha",
                  # Should be the eastern part of Rufiji
                  adm_div_2 == "Kibiti" ~ "Rufiji",
                  # created in 2016 from northen part of Bagamoyo 
                  adm_div_2 == "Chalinze" ~ "Bagamoyo",
                  # adm_div_1 Tanga
                  adm_div_2 == "Korogwe Rural" ~ "Korogwe",
                  # Created in 2013 from Lushoto
                  adm_div_2 == "Bumbuli" ~ "Lushoto",
                  # adm_div_1 Kilimanjaro
                  adm_div_2 == "Moshi Rural" ~ "Moshi",
                  .default = adm_div_2))

2.4 Crop the spatial variables

The spatial variables variables we have just load have a global coverage. It might be convenient to reduce the coverage to just the countries we are interested in. We can do this by using the terra::crop() function and the administrative divisions.

tpre_cntry <- purrr::map(adm_div, 
                         crop_with_buffer,
                         raster = tpre,
                         buffer = 1)

tevo_cntry <- purrr::map(adm_div, 
                         crop_with_buffer,
                         raster = tevo,
                         buffer = 1)

2.5 Weather Variable Transformation

# From meter to millimeters
pre_cntry_mm <- purrr::map(tpre_cntry, ~ .x*1000)

tevo_cntry_mm <- purrr::map(tevo_cntry, ~ .x*1000)

wb_cntry_mm <- purrr::map2(pre_cntry_mm, tevo_cntry_mm, \(x, y) x + y)

2.6 Extract

For the extraction of the weather variables, we use the function extract_cell_by_poly(). This function doesn’t aggregate the values within the polygons but extract all the cell values within the polygon separately. This is useful for us as we want to compute the SPI and SPEI for each cell and only later aggregate at the polygon level.

tpre_adm <- purrr::map2(pre_cntry_mm,
                        adm_div,
                        extract_cell_by_poly)

wb_adm <- purrr::map2(wb_cntry_mm,
                      adm_div,
                      extract_cell_by_poly)

2.7 Compute SPI and SPEI

We now compute the SPI and the SPEI with the function compute_spei() and compute_spi(). These functions requires the water balance and the precipitation time series for each location and the time scale at which the indices are computed.

To compute the SPEI, it is recommended to use at least 30 years of observation to ensure a good estimation of the parameters. More years can strength the estimation but the results can be affected by climate change: if there have been a change in the climate parameters, old observations might be not indicative of the current situation affecting the estimation. There are no clear rule on this, so we leave add the possibility to select the time range of observation with the function select_by_dates(). The function requires both or just one between the starting date, from, and the end date to. If both are provide the the function select between the two dates, if only from is provided the function selects all date after, and if only to is provided the function selects all date before.

Looking at the result, we see first is the ID column, that we will use to merge back with the survey. The other columns contain the SPEI observations over time specific to each coordinate.

spei6 <- purrr::map2(wb_adm,
                     names(wb_adm),
                     compute_spei,
                     time_scale = 6) |>
  purrr::map(agg_to_adm_div,
             match_col = "^X[0-9]")
Computing SPEI: BurkinaFaso22 
Computing SPEI: Ethiopia19 
Computing SPEI: Malawi19 
Computing SPEI: Mali22 
Computing SPEI: Nigeria18 
Computing SPEI: Tanzania21 
Computing SPEI: Togo22 
spei3 <- purrr::map2(wb_adm,
                     names(wb_adm),
                     compute_spei,
                     time_scale = 3) |>
  purrr::map(agg_to_adm_div,
             match_col = "^X[0-9]")
Computing SPEI: BurkinaFaso22 
Computing SPEI: Ethiopia19 
Computing SPEI: Malawi19 
Computing SPEI: Mali22 
Computing SPEI: Nigeria18 
Computing SPEI: Tanzania21 
Computing SPEI: Togo22 
spei1 <- purrr::map2(wb_adm,
                     names(wb_adm),
                     compute_spei,
                     time_scale = 1) |>
  purrr::map(agg_to_adm_div,
             match_col = "^X[0-9]")
Computing SPEI: BurkinaFaso22 
Computing SPEI: Ethiopia19 
Computing SPEI: Malawi19 
Computing SPEI: Mali22 
Computing SPEI: Nigeria18 
Computing SPEI: Tanzania21 
Computing SPEI: Togo22 
spi6 <- purrr::map2(tpre_adm,
                    names(tpre_adm),
                    compute_spi,
                    time_scale = 6) |>
  purrr::map(agg_to_adm_div,
             match_col = "^X[0-9]")
Computing SPI: BurkinaFaso22 
Computing SPI: Ethiopia19 
Computing SPI: Malawi19 
Computing SPI: Mali22 
Computing SPI: Nigeria18 
Computing SPI: Tanzania21 
Computing SPI: Togo22 
spi3 <- purrr::map2(tpre_adm,
                    names(tpre_adm),
                    compute_spi,
                    time_scale = 3) |>
  purrr::map(agg_to_adm_div,
             match_col = "^X[0-9]")
Computing SPI: BurkinaFaso22 
Computing SPI: Ethiopia19 
Computing SPI: Malawi19 
Computing SPI: Mali22 
Computing SPI: Nigeria18 
Computing SPI: Tanzania21 
Computing SPI: Togo22 
spi1 <- purrr::map2(tpre_adm,
                    names(tpre_adm),
                    compute_spi,
                    time_scale = 1) |>
  purrr::map(agg_to_adm_div,
             match_col = "^X[0-9]")
Computing SPI: BurkinaFaso22 
Computing SPI: Ethiopia19 
Computing SPI: Malawi19 
Computing SPI: Mali22 
Computing SPI: Nigeria18 
Computing SPI: Tanzania21 
Computing SPI: Togo22 

We have computed the standardized indicators for each cells but we need additional information to assign each cells to the administrative divisions. We take advantage of the column ID_adm_div to merge the standardized indices with the administrative divisions.

spei1_adm <- purrr::map(adm_div, terra::values) |>
  purrr::map2(spei1, merge_by_common)

spei3_adm <- purrr::map(adm_div, terra::values) |>
  purrr::map2(spei3, merge_by_common)

spei6_adm <- purrr::map(adm_div, terra::values) |>
  purrr::map2(spei6, merge_by_common)
spi1_adm <- purrr::map(adm_div, terra::values) |>
  purrr::map2(spi1, merge_by_common)

spi3_adm <- purrr::map(adm_div, terra::values) |>
  purrr::map2(spi3, merge_by_common)

spi6_adm <- purrr::map(adm_div, terra::values) |>
  purrr::map2(spi6, merge_by_common)

2.8 Merge with Survey

srvy_all_geo <- c(srvy_adm_geo, srvy_coord_geo)
srvy_all_geo <- srvy_all_geo[sort(names(srvy_all_geo))]

spi1_hh <- purrr::map2(srvy_all_geo,
                       spi1_adm,
                       merge_by_common) |>
    purrr::map(select_by_interview,
               interview = interview_date,
               interval = "2 year",
               wide = FALSE)
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!
spi3_hh <- purrr::map2(srvy_all_geo,
                       spi3_adm,
                       merge_by_common) |>
    purrr::map(select_by_interview,
               interview = interview_date,
               interval = "2 year",
               wide = FALSE)
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!
spi6_hh <- purrr::map2(srvy_all_geo,
                       spi6_adm,
                       merge_by_common) |>
    purrr::map(select_by_interview,
               interview = interview_date,
               interval = "2 year",
               wide = FALSE)
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!
spei1_hh <- purrr::map2(srvy_all_geo,
                        spei1_adm,
                        merge_by_common) |>
    purrr::map(select_by_interview,
               interview = interview_date,
               interval = "2 year",
               wide = FALSE)
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!
spei3_hh <- purrr::map2(srvy_all_geo,
                        spei3_adm,
                        merge_by_common) |>
    purrr::map(select_by_interview,
               interview = interview_date,
               interval = "2 year",
               wide = FALSE)
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!
spei6_hh <- purrr::map2(srvy_all_geo,
                        spei6_adm,
                        merge_by_common) |>
    purrr::map(select_by_interview,
               interview = interview_date,
               interval = "2 year",
               wide = FALSE)
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!
Missing interview are dropped!

2.9 Save

The final step of the code is to save the result. In this case, we save it as a dta file using the haven::write_dta() function.

spei6_hh |>
    dplyr::bind_rows() |>
    haven::write_dta(file.path(path_to_result, "spei6.dta"))
spei3_hh |>
    dplyr::bind_rows() |>
    haven::write_dta(file.path(path_to_result, "spei3.dta"))
spei1_hh |>
    dplyr::bind_rows() |>
    haven::write_dta(file.path(path_to_result, "spei1.dta"))


spi6_hh |>
    dplyr::bind_rows() |>
    haven::write_dta(file.path(path_to_result, "spi6.dta"))
spi3_hh |>
    dplyr::bind_rows() |>
    haven::write_dta(file.path(path_to_result, "spi3.dta"))
spi1_hh |>
    dplyr::bind_rows() |>
    haven::write_dta(file.path(path_to_result, "spi1.dta"))